
//
// Created by 王岩松 on 2020/3/31.
//

#include "pitch.h"

void Pitch::Pitch_create(double _physicalDuration, int _nFrames, double _dt, double _t1) {
    physicalDuration = _physicalDuration;
    nFrames = _nFrames;
    dt = _dt;
    t1 =_t1;
    candidates.resize(_nFrames);
}
int Pitch::Sampled_xToLowIndex(double time) {
    return floor((time - t1)/dt);
}
double Pitch::Sampled_indexToX(int index) {
    return t1 + dt * index;
}
double Pitch::getPitchAtX(double time, bool interpolate) {
    const double undefined = 0;

    if(time < 0 || time > physicalDuration) return undefined;
    if(interpolate) {
        const double index_real = (time - t1) / dt;
        const int leftIndex = floor(index_real);
        int nearIndex, farIndex;
        double phase = index_real - leftIndex;
        if(phase < 0.5) {
            nearIndex = leftIndex;
            farIndex = leftIndex + 1;
        } else {
            farIndex = leftIndex;
            nearIndex = leftIndex + 1;
            phase = 1.0 - phase;
        }
        if(nearIndex < 0 || nearIndex >= nFrames)
            return undefined;

        const double nearValue = candidates[nearIndex].frequnce[0];
        if(farIndex < 0 || farIndex >= nFrames)
            return nearValue;

        const double farValue = candidates[farIndex].frequnce[0];
        return nearValue + phase * (farValue - nearValue);
    }

    const double index_real = (time - t1) / dt;
    const int leftIndex = floor(index_real);
    int nearIndex;
    double phase = index_real - leftIndex;
    if(phase < 0.5) {
        nearIndex = leftIndex;
    } else {
        nearIndex = leftIndex + 1;
    }
    if(nearIndex < 0 || nearIndex >= nFrames)
        return undefined;

    const double nearValue = candidates[nearIndex].frequnce[0];
    return nearValue;
}
void Pitch::Free() {
    for(int i = 0; i < candidates.size(); ++ i) {
        vector<double>().swap(candidates[i].frequnce);
        vector<double>().swap(candidates[i].strength);
    }
    vector<Pitch_frame>().swap(candidates);
}


double NUM_interpolate_sinc (vector<double>& y, double x, int maxDepth) {
    const double PI = 3.1415926535897932384626433832795;
    int ix, midleft = (int) floor (x), midright = midleft + 1, left, right;
    double result = 0.0, a, halfsina, aa, daa;
    if (y.size() < 1) return 0;
    if (x > y.size()) return y [y.size()-1];
    if (x < 1) return y [1-1];
    if (x == midleft) return y [midleft-1];
    /* 1 < x < y.size && x not integer: interpolate. */
    if (maxDepth > midright - 1) maxDepth = midright - 1;
    if (maxDepth > y.size() - midleft) maxDepth = y.size() - midleft;
    if (maxDepth <= NUM_VALUE_INTERPOLATE_NEAREST) return y [(int) floor (x + 0.5)];
    if (maxDepth == NUM_VALUE_INTERPOLATE_LINEAR) return y [midleft-1] + (x - midleft) * (y [midright-1] - y [midleft-1]);
    if (maxDepth == NUM_VALUE_INTERPOLATE_CUBIC) {
        double yl = y [midleft], yr = y [midright];
        double dyl = 0.5 * (yr - y [midleft - 1 -1]), dyr = 0.5 * (y [midright + 1 -1] - yl);
        double fil = x - midleft, fir = midright - x;
        return yl * fir + yr * fil - fil * fir * (0.5 * (dyr - dyl) + (fil - 0.5) * (dyl + dyr - 2 * (yr - yl)));
    }
    left = midright - maxDepth;
    right = midleft + maxDepth;
    a = PI * (x - midleft);
    halfsina = 0.5 * sin (a);
    aa = a / (x - left + 1);
    daa = PI / (x - left + 1);
    for (ix = midleft; ix >= left; ix --) {
        double d = halfsina / a * (1.0 + cos (aa));
        result += y [ix -1] * d;
        a += PI;
        aa += daa;
        halfsina = - halfsina;
    }
    a = PI * (midright - x);
    halfsina = 0.5 * sin (a);
    aa = a / (right - x + 1);
    daa = PI / (right - x + 1); \
    for (ix = midright; ix <= right; ix ++) {
        double d = halfsina / a * (1.0 + cos (aa));
        result += y [ix -1] * d;
        a += PI;
        aa += daa;
        halfsina = - halfsina;
    }
    return result;
}





double NUMminimize_brent (double (*f) (double x, void *closure), double a, double b, void *closure, double tol, double *fx) {
    constexpr double NUM_goldenSection = 0.6180339887498948482045868343656381177203;
    double x, v, fv, w, fw;
    const double golden = 1 - NUM_goldenSection;
    const double sqrt_epsilon = sqrt (2.2204460492503131e-16);
    const int itermax = 60;

    assert (tol > 0 && a < b);

    /*
        First step - golden section
    */
    v = a + golden * (b - a);
    fv = (*f) (v, closure);
    x = v;
    w = v;
    *fx = fv;
    fw = fv;

    for (int iter = 1; iter <= itermax; iter ++) {
        const double middle_range = (a + b) / 2.0;
        const double tol_act = sqrt_epsilon * fabs (x) + tol / 3.0;
        double range = b - a;
        if (fabs (x - middle_range) + range / 2.0 <= 2.0 * tol_act)
            return x;

        // Obtain the golden section step

        double new_step = golden * (x < middle_range ? b - x : a - x);

        // Decide if the parabolic interpolation can be tried

        if (fabs (x - w) >= tol_act) {
            /*
                Interpolation step is calculated as p/q;
                division operation is delayed until last moment.
            */

            const double t = (x - w) * (*fx - fv);
            double q = (x - v) * (*fx - fw);
            double p = (x - v) * q - (x - w) * t;
            q = 2.0 * (q - t);

            if (q > 0.0)
                p = -p;
            else
                q = -q;

            /*
                If x+p/q falls in [a,b], not too close to a and b,
                and isn't too large, it is accepted.
                If p/q is too large then the golden section procedure can
                reduce [a,b] range.
            */

            if (fabs (p) < fabs (new_step * q) &&
                p > q * (a - x + 2.0 * tol_act) &&
                p < q * (b - x - 2.0 * tol_act)) {
                new_step = p / q;
            }
        }

        // Adjust the step to be not less than tolerance.

        if (fabs (new_step) < tol_act)
            new_step = new_step > 0.0 ? tol_act : - tol_act;

        // Obtain the next approximation to min	and reduce the enveloping range

        {
            const double t = x + new_step;	// Tentative point for the min
            const double ft = (*f) (t, closure);
            /*
                If t is a better approximation, reduce the range so that
                t would fall within it. If x remains the best, reduce the range
                so that x falls within it.
            */
            if (ft <= *fx) {
                if (t < x)
                    b = x;
                else
                    a = x;
                v = w;
                w = x;
                x = t;
                fv = fw;
                fw = *fx;
                *fx = ft;
            } else {
                if (t < x)
                    a = t;
                else
                    b = t;

                if (ft <= fw || w == x) {
                    v = w;
                    w = t;
                    fv = fw;
                    fw = ft;
                } else if (ft <= fv || v == x || v == w) {
                    v = t;
                    fv = ft;
                }
            }
        }
    }
    return x;
}
struct improve_params {
    int depth;
    vector<double> y;
    int isMaximum;
};
static double improve_evaluate (double x, void *closure) {
    struct improve_params *me = (struct improve_params *) closure;
    double y = NUM_interpolate_sinc (me->y, x, me-> depth);
    return me->isMaximum ? - y : y;
}
double NUMimproveExtremum (vector<double>& y, int ixmid, int interpolation, double *ixmid_real, bool isMaximum) {
    struct improve_params params;
    double result;
    if (ixmid <= 1) { *ixmid_real = 1; return y [1 -1]; }
    if (ixmid >= y.size()) { *ixmid_real = y.size(); return y [y.size()-1]; }
    if (interpolation <= NUM_PEAK_INTERPOLATE_NONE) { *ixmid_real = ixmid; return y [ixmid -1]; }
    if (interpolation == NUM_PEAK_INTERPOLATE_PARABOLIC) {
        double dy = 0.5 * (y [ixmid + 1 -1] - y [ixmid - 1 -1]);
        double d2y = 2 * y [ixmid -1] - y [ixmid - 1 -1] - y [ixmid + 1 -1];
        *ixmid_real = ixmid + dy / d2y;
        return y [ixmid -1] + 0.5 * dy * dy / d2y;
    }
    /* Sinc interpolation. */
    params. depth = interpolation == NUM_PEAK_INTERPOLATE_SINC70 ? 70 : 700;
    params. y = y;
    params. isMaximum = isMaximum;
    /*return isMaximum ?
        - NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate, & params, 1e-10, 1e-11, ixmid_real) :
        NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate, & params, 1e-10, 1e-11, ixmid_real);*/
    *ixmid_real = NUMminimize_brent (improve_evaluate, ixmid - 1, ixmid + 1, & params, 1e-10, & result);
    return isMaximum ? - result : result;
}

double NUMimproveMaximum (vector<double>& y, int ixmid, int interpolation, double *ixmid_real)
{ return NUMimproveExtremum (y, ixmid, interpolation, ixmid_real, true); }


int Pitch::Pitch_getMaxnCandidates() {
    int ans = 0;
    for(int i = 0; i < nFrames; ++ i) {
        ans = max(ans, candidates[i].nCandidates);
    }
    return ans;
}

inline static bool Pitch_util_frequencyIsVoiced (double f, double ceiling) {
    return f > 0.0 && f < ceiling;   // note: return false is f is NaN
}


void Pitch_pathFinder (Pitch &me, double silenceThreshold, double voicingThreshold,
                       double octaveCost, double octaveJumpCost, double voicedUnvoicedCost,
                       double ceiling, int pullFormants) {

    cout << "Pitch_pathFinder start" << endl;

    int nx = me.nFrames;
    double dx = me.dt;

    const int maxnCandidates = me.Pitch_getMaxnCandidates();
    cout << "maxnCandidates = " << maxnCandidates << endl;

    int place;
    volatile double maximum, value;
    const double ceiling2 = ( pullFormants ? 2.0 * ceiling : ceiling );
    /* Next three lines 20011015 */
    const double timeStepCorrection = 0.01 / dx;
    octaveJumpCost *= timeStepCorrection;
    voicedUnvoicedCost *= timeStepCorrection;

    vector<vector<double>> delta(nx, vector<double>(maxnCandidates,0));
    vector<vector<int>> psi(nx, vector<int>(maxnCandidates,0));

    cout << "delta start"<< endl;

    for (int iframe = 0; iframe < nx; iframe ++) {
        const Pitch_frame frame = me.candidates[iframe];
        double unvoicedStrength = ( silenceThreshold <= 0 ? 0.0 :
                                    2.0 - frame.intensity / (silenceThreshold / (1.0 + voicingThreshold)) );
        unvoicedStrength = voicingThreshold + std::max (0.0, unvoicedStrength);
        for (int icand = 0; icand < frame.nCandidates; icand ++) {

            const bool voiceless = ! Pitch_util_frequencyIsVoiced (frame.frequnce[icand], ceiling2);
            delta [iframe] [icand] = ( voiceless ? unvoicedStrength :
                                       frame.strength[icand] - octaveCost * log2 (ceiling / frame.frequnce[icand]) );
        }
    }

    cout << "delta end"<< endl;

    /* Look for the most probable path through the maxima. */
    /* There is a cost for the voiced/unvoiced transition, */
    /* and a cost for a frequency jump. */
    cout << "curDelta start"<< endl;

    for (int iframe = 1; iframe < nx; iframe ++) {
        const Pitch_frame &prevFrame = me.candidates[iframe - 1], &curFrame = me.candidates[iframe];
        const vector<double> &prevDelta = delta [iframe - 1];
        vector<double> &curDelta = delta [iframe];
        vector<int> &curPsi = psi [iframe];
        for (int icand2 = 0; icand2 < curFrame.nCandidates; icand2 ++) {
            const double f2 = curFrame.frequnce[icand2];
            maximum = -1e30;
            place = -1;
            for (int icand1 = 0; icand1 < prevFrame.nCandidates; icand1 ++) {
                double f1 = prevFrame.frequnce[icand1];
                double transitionCost;
                const bool previousVoiceless = ! Pitch_util_frequencyIsVoiced (f1, ceiling2);
                const bool currentVoiceless = ! Pitch_util_frequencyIsVoiced (f2, ceiling2);
                if (currentVoiceless) {
                    if (previousVoiceless) {
                        transitionCost = 0.0;   // both voiceless
                    } else {
                        transitionCost = voicedUnvoicedCost;   // voiced-to-unvoiced transition
                    }
                } else {
                    if (previousVoiceless) {
                        transitionCost = voicedUnvoicedCost;   // unvoiced-to-voiced transition
                    } else {
                        transitionCost = octaveJumpCost * fabs (log2 (f1 / f2));   // both voiced
                    }
                }
                value = prevDelta [icand1] - transitionCost + curDelta [icand2];
                //if (Melder_debug == 33) Melder_casual (U"Frame ", iframe, U", current candidate ", icand2,
                //  U" (delta ", curDelta [icand2], U"), previous candidate ", icand1, U" (delta ", prevDelta [icand1], U"), ",
                //	U"transition cost ", transitionCost, U", value ", value, U", maximum ", maximum);
                if (value > maximum) {
                    maximum = value;
                    place = icand1;
                }
            }
            curDelta [icand2] = maximum;
            curPsi [icand2] = place;
        }
    }
    cout << "curDelta end"<< endl;

    /* Find the end of the most probable path. */

    place = 0;
    maximum = delta [nx-1] [place];
    for (int icand = 1; icand < me.candidates[nx-1].nCandidates; icand ++) {
        if (delta [nx-1] [icand] > maximum) {
            place = icand;
            maximum = delta [nx-1] [place];
        }
    }
    cout << "find place end"<< endl;

    /* Backtracking: follow the path backwards. */

    for (int iframe = nx-1; iframe >= 0; iframe --) {
        Pitch_frame &frame = me.candidates[iframe];
        std::swap (frame.frequnce[0], frame.frequnce[place]);
        std::swap (frame.strength[0], frame.strength[place]);
        place = psi [iframe] [place];
//        cout << place << " < ";
    }
    cout << "swap place end"<< endl;

    /* Pull formants: devoice frames with frequencies between ceiling and ceiling2. */

    if (ceiling2 > ceiling) {
        for (int iframe = nx-1; iframe >= 0; iframe --) {
            Pitch_frame &frame = me.candidates[iframe];
            double & winner_feq = frame.frequnce[0];
            double & winner_str = frame.strength[0];

            const double f = winner_feq;
            if (f > ceiling && f < ceiling2) {
                for (int icand = 1; icand < frame.nCandidates; icand ++) {
                    double &loser_feq = frame.frequnce[icand];
                    double &loser_str = frame.strength[icand];
                    if (loser_feq == 0.0) {
                        std::swap (winner_feq, loser_feq);
                        std::swap (winner_str, loser_str);
                        break;
                    }
                }
            }
        }
    }
    Pitch::freeVV(delta);
    Pitch::freeVV(psi);
}

static void Sound_into_PitchFrame (sound me, Pitch_frame &pitchFrame, double t,
                                   double minimumPitch, int maxnCandidates, int method, double voicingThreshold, double octaveCost,
                                    double dt_window, int nsamp_window, int halfnsamp_window,
                                   int maximumLag, int nsampFFT, int nsamp_period, int halfnsamp_period,
                                   int brent_ixmax, int brent_depth, double globalPeak,
                                   vector<double> & frame, vector<double> & ac, vector<double> & window, vector<double> & windowR,
                                   double *r, vector<int> & imax, double & localMean){

    int nx = me.getlength();
    double dx = 1.0/me.getfrequency();

    int leftSample = me.Sampled_xToLowIndex (t), rightSample = leftSample + 1;
    int startSample, endSample;

    startSample = rightSample - nsamp_period;
    endSample = leftSample + nsamp_period;
    assert (startSample >= 0);
    assert (endSample < nx);
    localMean = 0.0;
    for (int i = startSample; i <= endSample; i ++)
        localMean += me[i];
    localMean /= 2 * nsamp_period;

    /*
        Copy a window to a frame and subtract the local mean.
        We are going to kill the DC component before windowing.
    */
    startSample = rightSample - halfnsamp_window;
    endSample = leftSample + halfnsamp_window;
    assert (startSample >= 0);
    assert (endSample <= nx);

    for (int j = 0, i = startSample; j < nsamp_window; j ++)
        frame [j] = (me[i ++] - localMean ) * window [j];
    for (int j = nsamp_window; j < nsampFFT; j ++)
        frame [j] = 0.0;

    /*
        Compute the local peak; look half a longest period to both sides.
    */
    double localPeak = 0.0;
    if ((startSample = halfnsamp_window + 1 - halfnsamp_period) < 1)
        startSample = 0;
    if ((endSample = halfnsamp_window + halfnsamp_period) > nsamp_window)
        endSample = nsamp_window-1;

    for (int j = startSample; j <= endSample; j ++) {
        double value = fabs (frame[j]);
        if (value > localPeak)
            localPeak = value;
    }

    int & nCandidates = pitchFrame.nCandidates;
    nCandidates = 0;
    double & intensity = pitchFrame.intensity;
    intensity = ( localPeak > globalPeak ? 1.0 : localPeak / globalPeak );

/*
            The FFT of the autocorrelation is the power spectrum.
        */
    for (int i = 0; i < nsampFFT; i ++)
        ac [i] = 0.0;

    init_rfft(nsampFFT);
    rfft(frame,nsampFFT);   // complex spectrum
    ac [0] += frame[0] * frame[0];   // DC component
    for (int i = 1; i < nsampFFT-1; i += 2)
        ac [i] += frame [i] * frame [i] + frame [i+1] * frame [i+1];   // power spectrum
    ac [nsampFFT-1] += frame [nsampFFT-1] * frame [nsampFFT-1];   // Nyquist frequency

    irfft(ac,nsampFFT);   // autocorrelation

    /*
        Normalize the autocorrelation to the value with zero lag,
        and divide it by the normalized autocorrelation of the window.
    */
    r [0] = 1.0;
    for (int i = 1; i <= brent_ixmax; i ++)
        r [- i] = r [i] = ac [i] / (ac [0] * windowR [i]);

    pitchFrame.frequnce.resize (nCandidates = 1);
    pitchFrame.strength.resize (nCandidates = 1);
    pitchFrame.frequnce[0] = pitchFrame.strength[0] = 0.0;

    if (localPeak == 0.0)
        return;

    auto constVEC = [](double* a, int n) {
        vector<double> res(n);
        for(int i = 0; i < n; ++ i) {
            res[i] = a[i];
        }
        return res;
    };

    imax [1] = 0;
    for (int i = 2; i < maximumLag && i < brent_ixmax; i ++)
        if (r [i] > 0.5 * voicingThreshold &&   // not too unvoiced?
            r [i] > r [i-1] && r [i] >= r [i+1])   // maximum?
        {
            int place = 0;
            const double dr = 0.5 * (r [i+1] - r [i-1]), d2r = 2.0 * r [i] - r [i-1] - r [i+1];
            const double frequencyOfMaximum = 1.0 / dx / (i + dr / d2r);
            const int offset = - brent_ixmax - 1;
            vector<double> tt = constVEC(& r [offset + 1], brent_ixmax - offset);
            double strengthOfMaximum = NUM_interpolate_sinc (tt, 1.0 / dx / frequencyOfMaximum - offset, 30);

            if (strengthOfMaximum > 1.0)
                strengthOfMaximum = 1.0 / strengthOfMaximum;

            if (nCandidates < maxnCandidates) {   // is there still a free place?
                pitchFrame.frequnce. resize (++ nCandidates);
                pitchFrame.strength. resize (nCandidates);
                place = nCandidates;
            } else {
                /*
                    Try the place of the weakest candidate so far.
                */
                double weakest = 2.0;
                for (int iweak = 2; iweak <= maxnCandidates; iweak ++) {
                    /*
                        High frequencies are to be favoured
                        if we want to analyze a perfectly periodic signal correctly.
                    */
                    double localStrength = pitchFrame. strength[iweak -1] - octaveCost *log2 (minimumPitch / pitchFrame.frequnce[iweak -1]);
                    if (localStrength < weakest) {
                        weakest = localStrength;
                        place = iweak;
                    }
                }
                /*
                    If this maximum is weaker than the weakest candidate so far, give it no place.
                */
                if (strengthOfMaximum - octaveCost * log2 (minimumPitch / frequencyOfMaximum) <= weakest)
                    place = 0;
            }
            if (place) {   // have we found a place for this candidate?
                pitchFrame.frequnce[place -1] = frequencyOfMaximum;
                pitchFrame.strength[place -1] = strengthOfMaximum;
                imax [place] = i;
            }

        }

    /*
    Second pass: for extra precision, maximize sin(x)/x interpolation ('sinc').
*/
    vector<double> ttt;
    for (int i = 2; i <= nCandidates; i ++) {
        if (method != AC_HANNING || pitchFrame.frequnce[i -1] > 0.0 / dx) {
            double xmid, ymid;
            const int offset = - brent_ixmax - 1;
            ttt = constVEC (& r [offset + 1], brent_ixmax - offset);
            ymid = NUMimproveMaximum (ttt, imax [i] - offset,
                                      pitchFrame.frequnce[i -1] > 0.3 / dx ? NUM_PEAK_INTERPOLATE_SINC700 : brent_depth, & xmid);
            xmid += offset;
            pitchFrame.frequnce[i -1] = 1.0 / dx / xmid;
            if (ymid > 1.0)
                ymid = 1.0 / ymid;
            pitchFrame.strength [i -1] = ymid;
        }
    }
    Pitch::freeVector(ttt);
}

Pitch Sound_to_Pitch_any (sound me,
                          double dt, double minimumPitch, double periodsPerWindow, int maxnCandidates,
                          int method,
                          double silenceThreshold, double voicingThreshold,
                          double octaveCost, double octaveJumpCost, double voicedUnvoicedCost, double ceiling){
    const double PI = acos(-1.0);
    double t1;
    int numberOfFrames;
    int nsampFFT;
    double interpolation_depth;
    int brent_ixmax, brent_depth;
    double globalPeak;

    try {
        if (maxnCandidates < ceiling / minimumPitch)
            maxnCandidates = floor (ceiling / minimumPitch);
        if (dt <= 0.0)
            dt = periodsPerWindow / minimumPitch / 4.0;

//        switch (method) {
//            case AC_HANNING:
                brent_depth = NUM_PEAK_INTERPOLATE_SINC70;
                interpolation_depth = 0.5;
//                break;
//        }

        double Fs = me.getfrequency();
        double dx = 1.0 / me.getfrequency();
        int nx = me.getlength();

        double duration = 1.0/dx * nx;
        if (minimumPitch < periodsPerWindow / duration)
            throw ("To analyse this Sound, “minimum pitch” must not be less than "+to_string(periodsPerWindow / duration)+" Hz.");
        int nsamp_period = floor (me.getfrequency() / minimumPitch);
        int halfnsamp_period = nsamp_period / 2 + 1;
        if(ceiling > 0.5 * me.getfrequency())
            ceiling = 0.5 * me.getfrequency();

        double dt_window = periodsPerWindow / minimumPitch;
        int nsamp_window = floor (dt_window * me.getfrequency());
        int halfnsamp_window = nsamp_window / 2 - 1;
        if (halfnsamp_window < 2)
            throw ("Analysis window too short.");
        nsamp_window = halfnsamp_window * 2;
        const int minimumLag = std::max (2.0, floor (me.getfrequency() / ceiling));
        const int maximumLag = std::min (floor (nsamp_window / periodsPerWindow) + 2, nsamp_window * 1.0);
        const double physicalDuration = 1.0 * nx / Fs;
        // nFrames , t1
        numberOfFrames = 1 + floor ((physicalDuration - dt_window) / dt);
        t1 = 0.0 + 0.5 * (physicalDuration - 1.0/Fs - (numberOfFrames - 1) * dt);   // centre of first frame

        Pitch thee; thee.Pitch_create(physicalDuration,numberOfFrames,dt,t1);
        for(int i = 0; i < numberOfFrames; ++ i) {
            thee.candidates[i].frequnce.resize(maxnCandidates+1);
            thee.candidates[i].strength.resize(maxnCandidates+1);
        }
        cout << "Pitch_create" << endl;

        globalPeak = 0.0;
        auto cal_mean = [&](sound &me) {
            int len = me.getlength();
            long double avg = 0;
            for(int i = 0; i < len; ++ i)
                avg += 1.0 * me[i];
            return (double)(avg / len);
        };

        const double mean = cal_mean(me);
        for (int i = 0; i < nx; i ++) {
                double value = fabs (me[i] - mean);
                if (value > globalPeak)
                    globalPeak = value;
        }

        cout << "globalPeak = " << globalPeak << endl;

        if (globalPeak == 0.0)
            return thee;

        vector<double> window, windowR;

        nsampFFT = 1;
        while (nsampFFT < nsamp_window * (1 + interpolation_depth))
            nsampFFT *= 2;

        windowR. resize (nsampFFT);
        window. resize (nsamp_window);

        for (int i = 1; i <= nsamp_window; i ++)
            window [i-1] = 0.5 - 0.5 * cos (i * 2 * PI / (nsamp_window + 1));
        for (int i = 0; i < nsamp_window; i ++)
            windowR [i] = window [i];

        cout << "window build :" << endl;
//        for(auto aa: window) cout << aa << " "; cout << "\n";

        init_rfft(nsampFFT);
        cout << "init_fft" << endl;
        rfft(windowR, nsampFFT);
        cout << "fft" << endl;

        windowR [0] *= windowR [0];   // DC component
        for (int i = 1; i < nsampFFT-1; i += 2) {
            windowR [i] = windowR [i] * windowR [i] + windowR [i + 1] * windowR [i + 1];
            windowR [i + 1] = 0.0;   // power spectrum: square and zero
        }
        windowR [nsampFFT-1] *= windowR [nsampFFT-1];   // Nyquist frequency

        cout << "ifft_start" << endl;
        irfft(windowR, nsampFFT);
        cout << "ifft_end" << endl;

        for (int i = 1; i <= nsamp_window-1; i ++)
            windowR [i] /= windowR [0];   // normalize
        windowR [0] = 1.0;   // normalize

        cout << "windowR build: " << endl;
//        for(auto aa: windowR) cout << aa << " "; cout << "\n";

        brent_ixmax = floor (nsamp_window * interpolation_depth);
        cout << "brent_ixmax = " << brent_ixmax << endl;

        vector<double> frame(nsampFFT), ac(nsampFFT);

        double * buffer = new double[2 * nsamp_window + 1];
        double * r = &buffer[1 + nsamp_window];
        vector<int> imax(maxnCandidates+1);
        double localMean;

        for (int iframe = 0; iframe < numberOfFrames; iframe ++) {
            const double t = thee.Sampled_indexToX(iframe);
            Sound_into_PitchFrame (me, thee.candidates[iframe], t, minimumPitch, maxnCandidates,
                    method,voicingThreshold,octaveCost,dt_window,nsamp_window,halfnsamp_window,
                    maximumLag,nsampFFT,nsamp_period,halfnsamp_period,brent_ixmax,brent_depth,globalPeak,
                    frame, ac, window, windowR, r, imax, localMean);
        }

        delete(buffer);

        Pitch_pathFinder (thee, silenceThreshold, voicingThreshold,
                          octaveCost, octaveJumpCost, voicedUnvoicedCost, ceiling, false);

        Pitch::freeVector(imax);
        Pitch::freeVector(frame);
        Pitch::freeVector(ac);
        Pitch::freeVector(window);
        Pitch::freeVector(windowR);

        return thee;

    }catch (string s) {
        cerr << s << endl;
    }
}

Pitch Sound_to_Pitch_ac (sound me,
                         double dt, double minimumPitch, double periodsPerWindow, int maxnCandidates, int accurate,
                         double silenceThreshold, double voicingThreshold,
                         double octaveCost, double octaveJumpCost, double voicedUnvoicedCost, double ceiling)
{
    return Sound_to_Pitch_any (me, dt, minimumPitch, periodsPerWindow, maxnCandidates, accurate,
                               silenceThreshold, voicingThreshold, octaveCost, octaveJumpCost, voicedUnvoicedCost, ceiling);
}

Pitch Sound_to_Pitch (sound me, double timeStep, double minimumPitch, double maximumPitch) {
    return Sound_to_Pitch_ac (me, timeStep, minimumPitch,
                              3.0, 15, false, 0.03, 0.45, 0.01, 0.35, 0.14, maximumPitch);
}

